Multidimensional replica-exchange method for free-energy calculations 
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We have developed a new simulation algorithm for free-energy calculations. The method is a 
multidimensional extension of the replica-exchange method. While pairs of replicas with different 
j^ , temperatures are exchanged during the simulation in the original replica-exchange method, pairs 

of replicas with different temperatures and/or different parameters of the potential energy are ex- 
changed in the new algorithm. This greatly enhances the sampling of the conformational space and 
'■^ , allows accurate calculations of free energy in a wide temperature range from a single simulation run, 

^ ' using the weighted histogram analysis method. 

o 

_CJ . 

I. INTRODUCTION 

1—i '. 

J^ ] In complex systems such as a system of proteins, it is difficult to obtain accurate canonical distributions at low 

i^ ■ temperatures by the conventional molecular dynamics (MD) or Monte Carlo (MC) simulations. This is because 

I ' there exist a huge number of local-minimum states in the potential energy surface, and the simulations tend to get 

^\ . trapped in one of the local-minimum states. One popular way to overcome this difficulty is to perform a generalized- 

C^ ensemble simulation, which is based on non-Boltzmann probability weight factors so that a random walk in energy 

^D , space may be realized (for a review see |1[). The random walk allows the simulation to go over any energy barrier 

^p ■ and sample much wider configurational space than by conventional methods. Monitoring the energy in a single 

simulation run, one can obtain not only the global-minimum-energy state but also any thermodynamic quantities as 

a function of temperature for a wide temperature range. The latter is made possible by the single-histogram pi or 

multiple-histogram reweighting techniques (an extension of the multiple-histogram method is also referred to as 

1-^ , the weighted histogram analysis method (WHAM) Q). 

(— I ' Three of the most well-known generalized-ensemble methods are perhaps multicanonical algorithm [^, simulated 

Q tempering [pp], and replica-exchange method [pi- |12L (The replica-exchange method is also referred to as replica Monte 

O Garlo method ||9|, multiple Markov chain method |llj, and parallel tempering [[l2[.) These algorithms have already 

^ ^ been used in many applications in protein and related systems (see, for instance, Refs. p3| - p6[ for multicanonical 

•""j ■ algorithm, Refs. [^ - ^^1 foi' simulated tempering, and Refs. [Q - |^^ for replica-exchange method). 

KS The replica-exchange method (REM) has been drawing much attention recently because the probability weight 

^ , factors are essentially known a priori, whereas they are not in most other generalized-ensemble algorithms (and have 

■ ■ ■ ' to be determined by a tedius procedure). In REM the generalized ensemble consists of noninteracting copies (or 

replicas) of the original system with different temperatures. During a parallel MD or MC simulation of each replica. 
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a pair of replicas are exchanged every few steps. This procedure enforces random walks in the replica (temperature) 
space. 

In a previous work p4[ we worked out the details for the replica-exchange molecular dynamics algorithm. In this 
article we present a multidimensional extension of the replica-exchange method (similar generalizations of REM can 
also be found in Rcfs. M,E^). In the new algorithm, pairs of replicas with different temperatures and/or different 
parameters of the potential energy are exchanged. As an example of the applications of the multidimensional replica- 
exchange method, we discuss free-e nerg y calculations in detail. 

The umbrella sampling method |l8| and free-energy perturbation method, which is a special case of umbrella 
sampling, have been widely used to calculate the free energies in chemical processes |M - [M. In the umbrella 
sampling method, a reaction coordinate is chosen and free-energy profiles along the reaction coordinate are calculated. 
A series of independent simulations are performed to sample the relevant range of the coordinate. To cover the entire 
range of the coordinate, biasing potentials, which are called "umbrella potentials," are imposed. Thus, the system is 
restrained to remain near the prechosen value of the reaction coordinate specified by each umbrella potential, and a 
series of simulations with different umbrella potentials are performed. WHAM [Q is often employed to calculate the 
free-energy profiles from the histograms obtained by each simulation. 

Although the effectiveness of the umbrella sampling method is well known, its successful implementation requires 
a careful fine tuning. For instance, the choice of the umbrella potentials is very important. If the potentials are too 
strong, the conformational space sampled by each simulation becomes quite narrow. If the potentials are too weak, 
on the other hand, the system does not remain near the prechosen value of the reaction coordinate. The values of 
the coupling parameters A for the umbrella potentials should also be carefully chosen. Various generalizations of the 
umbrella sampling method have thus been introduced to sample the potential energy surface more effectively. The 
X-dynamics |5q] - |60[ is such an example, where the coupling parameter A is treated as a dynamical variable. Another 
example is the multicanonical WHAM pq], which combines the umbrella sampling with multicanonical algorithm. In 
the present article we develop yet another generalization of the umbrella sampling method (we refer to this method 
as replica- exchange umbrella sampling), which is based on the multidimensional extension of the replica-exchange 
method. 

In Section II the multidimensional extension of the replica-exchange method is described in detail. In particular, 
the replica-exchange umbrella sampling method is introduced. In Section III the results of the application of replica- 
exchange umbrella sampling to a blocked alanine trimer are given. Section IV is devoted to conclusions. 

II. METHODS 

Before we describe the multidimensional replica-exchange method (MREM), let us briefly review the original replica- 
exchange method (REM) g- [^ (see Ref. @ for details). 

We consider a system of N atoms with their coordinate vectors and momentum vectors denoted by 5 = {Qi, • • • , g^y} 
and p = {Pi, ■ ■ ■ jPat}, respectively. The Hamiltonian H(q,p) of the system is the sum of the kinetic energy K{p) and 
the potential energy E{q): 

H{q,p)^K{p)+Eiq) . (1) 

In the canonical ensemble at temperature T each state x = {q,p) with the Hamiltonian H(q,p) is weighted by the 
Boltzmann factor: 

Weix) = e-l^"^-^'P^ , (2) 

where the inverse temperature /3 is defined by /? = l/ksT (ks is Boltzmann's constant). 

The generalized ensemble for REM consists of M noninteracting copies (or, replicas) of the original system in 
the canonical ensemble at M different temperatures Tm {m = 1, • ■ • , M). We arrange the replicas so that there is 
always exactly one replica at each temperature. Then there is a one-to-one correspondence between replicas and 
temperatures; the label i (i = 1, ■ ■ ■ , M) for replicas is a permutation of the label m (m = 1, • ■ • , M) for temperatures, 
and vice versa: 



i = i{m) = f{m) , 
m = m{i) = f^^{i) 



(3) 



where f{m) is a permutation function of m and / ^{i) is its inverse. 



Let X = <^ Xi , • ■ • ,x^) f ~ 1 •^mi'i')' ■ ' ■ '•^mi'M') r stand for a "state" in this generalized ensemble. Here, the 

superscript and the subscript in xhl label the replica and the temperature, respectively. The state X is specified by 
the M sets of coordinates gt*' and momenta pW of N atoms in replica i at temperature T™ : 

xW^(gW,pH) . (4) 

Because the replicas are noninteracting, the weight factor for the state X in this generalized ensemble is given by the 
product of Boltzmann factors for each replica (or at each temperature) : 

Whem{X) = exp I -f]/3„(,oi? (g'^p''!) } = cxp I -f] /3™i/ (qW'")l,pW™)l) I , (5) 

where i{m) and m^i) are the permutation functions in Eq. (y). 

We now consider exchanging a pair of replicas in the generalized ensemble. Suppose we exchange replicas i and j 
which are at temperatures T™ and r„, respectively: 

The exchange of replicas can be written in more detail as 






(7) 



where the momenta are uniformly rescaled according to p4| 

V ^ m 

V -^n 



(8) 



In order for this exchange process to converge toward the equilibrium distribution based on Eq. (m, it is sufficient 
to impose the detailed balance condition on the transition probability w{X -^ X'): 

WREMiX) wiX ^ X') = Wrem{X') w{X' ^ X) . (9) 

From Eqs. (|), (|), (|), and (§), we have 

w(X ^ X') , ^, 

where 

A = /3„ (i? (g[^l) - i? (gW)) - /3„ (i? (g[^l) - E (gW)) , (11) 

= (/3™-/3„)(i?(-7[^l)-i?(gW)) . (12) 

This can be satisfied, for instance, by the usual Metropolis criterion pll : 

.(.Y^A-)..(.a|4i) = {Lp(-A),f:i;l!: m 

Note that because of the velocity rescaling of Eq. (g) the kinetic energy terms are cancelled out in Eqs. (|l^) (and 
(O)) and that the same criterion, Eqs. ( [12|) and (|l3|), which was originally derived for the Monte Carlo algorithm Q- 
| [12[ is recovered [ p4| . 

A simulation of the replica-exchange method (REM) |^- |12(| is then realized by alternately performing the following 
two steps: 



1. Each replica in canonical ensemble of the fixed temperature is simulated simultaneously and independently for 
a certain MC or MD steps. 

2. A pair of replicas, say Xm and x„ , are exchanged with the probability w (xm Xn ) in Eq. (Oh. 

In the present work, we employ the molecular dynamics algorithm for step 1. Note that in step 2 we exchange only 
pairs of replicas corresponding to neighboring temperatures, because the acceptance ratio of the exchange decreases 
exponentially with the difference of the two /3's (see Eqs. (|lj) and (|3|)). Note also that whenever a replica exchange 
is accepted in step 2, the permutation functions in Eq. (ra) are updated. 

The major advantage of REM over other generalized-ensemble methods such as multicanonical algorithm [^j and 
simulated tempering BM] lies in the fact that the weight factor is a priori known (see Eq. (0)), while in the latter 
algorithms the determination of the weight factors can be very tedius and time-consuming. A random walk in 
"temperature space" is realized for each replica, which in turn induces a random walk in potential energy space. This 
alleviates the problem of getting trapped in states of energy local minima. 

We now present our multidimensional extension of REM, which we refer to as the multidimensional replica- exchange 
method (MREM). The crucial observation that led to the new algorithm is: As long as we have M noninteracting 
replicas of the original system, the Hamiltonian H{q,p) of the system does not have to be identical among the replicas 
and it can depend on a parameter with different parameter values for different replicas. Namely, we can write the 
Hamiltonian for the i-th replica at temperature T,n as 

^m(g'^p'^l) - i^(/l) -i-i?A,„(g[*l) , (14) 

where the potential energy E\^^ depends on a parameter A™ and can be written as 

E^M^)^Eo{q^'^) + XmV{q^'^) . (15) 

This expression for the potential energy is often used in simulations. For instance, in umbrella sampling [J48| , £'o('Z) 
and V{q) can be respectively taken as the original potential energy and the "biasing" potential energy with the 
coupling parameter A,„. In simulations of spin systems, on the other hand, i?o(<z) and V{q) (here, q stands for spins) 
can be respectively considered as the zero-field term and the magnetization term coupled with the external field Am. 
While replica i and temperature Tm are in one-to-one correspondence in the original REM, replica i and "parameter 
set" Am = (Tm,Am) are in one-to-one correspondence in the new algorithm. Hence, the present algorithm can be 
considered as a multidimensional extension of the original replica-exchange method where the "parameter space" is 
one dimensional (i.e., Am = Tm)- Because the replicas are noninteracting, the weight factor for the state X in this 
new generalized ensemble is again given by the product of Boltzmann factors for each replica (see Eq. (||)): 

Wmrem(.X) = exp I -y0mU)HmU) (^I^p''') ^ = exp << - V p^Hm (qW'")],^^™)]) \ , (16) 



i=l 




where i{m) and m(i) are the permutation functions in Eq. (P). Then the same derivation that led to the original 
replica-exchange criterion (Eq. (|l3|)) follows, and we have the following transition probability of replica exchange (see 
E<t (0)): 



where 



A 



/3m [e^^ (gl-'l) - i?A„ (gW) ) - Pn {ex,^ (gl-'l) - i?A„ (g['l) ) • (18) 



Here, E\^ and E\^^ are the total potential energies (see Eq. (|15|)). Note that we need to newly evaluate the potential 
energy for exchanged coordinates, Ex^iq^-'^) and iJA„ (?'*'), because Ex^ and Ex„ are in general different functions. 
The method is particularly suitable for parallel computers. Because one can minimize the amount of information 
exchanged among nodes, it is best to assign each replica to each node (exchanging Tm, Ex^ and T„, Ex^ among 
nodes is much faster than exchanging coordinates and momenta). This means that we keep track of the permutation 
function m{i;t) = f~^{i;t) in Eq. (y) as a function of MD step t throughout the simulation. 

For obtaining the canonical distributions, the weighted histogram analysis method (WHAM) H] is particularly suit- 
able. Suppose we have made a single run of the present replica-exchange simulation with M replicas that correspond 



to M different parameter sets A^ = (Tm, Am) (w = 1, • ■ • , M). Let Nm{EQ,V) and n™ be, respectively, the potential- 
energy histogram and the total number of samples obtained for the m-th parameter set Am. The expectation value 
of a physical quantity A for any potential-energy parameter value A at any temperature T = l/fcs/3 is then given by 



Si 



where 



Y, MEo,V)PtAEo,V) 



< A >T.. 



Eo,V 



J2Pt,x{Eo,V) 

Eo,V 



(19) 



PT.xiEo,V) 



M 



Y,97nNAE^,V) 

771—1 



-PE^ 



(20) 



and 



,-/,r 



^Pt„,a„(^o,V^) 



(21) 



Ea,V 



Here, gm = 1 + 2r,„, and t,„ is the integrated autocorrelation time at temperature T„j with the parameter value A,„. 
Note that the unnormalized probability distribution Pt.x{Eq, V) and the "dimensionless" Helmholtz free energy /„ 
in Eqs. (gO|) and ( pi|) are solved self-consistently by iteration [^|j]. 

We can use this new replica-exchange method for free-energy calculations. We first describe the free-energy pertur- 
bation case. The potential energy is given by 



Ex{q) = Ei{q) + \{Ep{q) - Ei{q)) 



(22) 



where Ei and Ep are the potential energy for a "wild-type" molecule and a "mutated" molecule, respectively. Note 
that this equation has the same form as Eq. (p^). 

Our replica-exchange simulation is performed for M replicas with M different values of the parameters A„i = 
(Tm, Am). Since E\=o(,q) — Ei{q) and E\^i{q) = Epi^q), we should choose enough Am values distributed in the range 
between and 1 so that we may have sufficient replica exchanges. From the simulation, M histograms Nj^^Ej^ Ep — 
Ei), or equivalently Nm{Ei, Ep), are obtained. The Helmholtz free-energy difference of "mutation" at temperature 
T, AF = F\=i — Fx^o, can then be calculated from 



exp(-/3Ai^) 



Zt,\=o 



Y Pt.x=i{Ei,Ep) 

Ej ^E p 



(23) 



where Pt.\{Ej,Ep) are obtained from the WHAM equations of Eqs. (|20|) and (pi]). 

We now describe another free-energy calculation based on MREM applied to umbrella sampling jBq] , which we refer 
to as replica-exchange umbrella sampling (REUS). The potential energy is a generalization of Eq. (|l5| ) and is given by 



E^{q)^Eo{q)+J2^^'^VdQ) 



(24) 



where Eo{q) is the original unbiased potential, Ve{q) (i — 1, • • • ,L) are the biasing (umbrella) potentials, and A^^^ 
are the corresponding coupling constants (A — (A^^-', • ■ • , A'^^-')). Introducing a "reaction coordinate" ^, the umbrella 
potentials are usually written as harmonic restraints: 



VM)^Mi{<i)-M^ , (^ = i,---,L) 



(25) 



where di are the midpoints and ki are the strengths of the restraining potentials. We prepare M replicas with M 
different values of the parameters Am — (Tm,Am), and the replica-exchange simulation is performed. Since the 



umbrella potentials Ve {q) in Eq. (^5|) are all functions of the reaction coordinate ^ only, we can take the histogram 
Nrn{Eo, instead of Njn{Eo,Vi,- ■ ■ ,Vl). The WHAM equations of Eqs. ^) and (|l|) can then be written as 



^rA(^o,0 = 



M 



and 



J29;n'NmiEo,0 

2^ nmgm e '^" 






-PE^ 



(26) 



(27) 



The expectation value of a physical quantity A is now given by (see Eq. (09)) 



< A >T.X = 



Y,A{Eo.OPT,\iEo,0 

Ea£ 

Eo,i 



(28) 



The potential of mean force (PMF), or free energy as a function of the reaction coordinate, of the original, unbiased 
system at temperature T is given by 



^TXM^}iO = ~kBT\n 






(29) 



where {0} = (0, • • • , 0). In the examples presented below, replicas were chosen so that the potential energy for each 
replica includes exactly one (or zero) biasing potential. 



III. RESULTS AND DISCUSSION 



One of the applications of MREM, replica- exchange umbrella sampling (REUS), was tested for the system of a 
blocked peptide, alanine trimer. The N and C termini of the peptide were blocked with acetyl and N-methyl groups, 
respectively. Since the thermodynamic behavior of this peptide was extensively studied by the conventional umbrella 
sampling |^1|1 , it is a good test case to examine the effectiveness of the new method. All calculations were based on 
MD simulations, and the force field parameters were taken from the all-atom version of AMBER ^^ with a distance- 
dependent dielectric, e — r, which mimics the presence of solvent. The computer code developed in Refs. |p6| , |63| , 
which is based on Version 2 of PRESTO IQ, was used. The temperature during the MD simulations was controlled 
by the constraint method [£5[g6|. The unit time step was set to 0.5 fs, and we made an MD simulation of 4 x 10^ 
time steps (or 2.0 ns) for each replica, starting from an extended conformation. The data were stored every 20 steps 
(or 10 fs) for a total of 2 x 10^ snapshots. (Before taking the data, we made regular canonical MD simulations of 100 
ps for thermalization. For replica-exchange simulations an additional REM simulation of 100 ps was made for further 
thermalization. ) 

In Table I we summarize the parameters characterizing the replicas for the simulations performed in the present 
work. They are one original replica-exchange simulation (REMl), two replica-exchange umbrella sampling simulations 
(REUSl and REUS2), and two conventional umbrella sampling simulations (USl and US2). The purpose of the present 
simulations is to test the effectiveness of the replica-exchange umbrella sampling with respect to the conventional 
umbrella sampling (REUSl and REUS2 versus USl and US2). The original replica-exchange simulation without 
umbrella potentials (REMl) was also made to set a reference standard for comparison. For REMl, replica exchange 
was tried every 20 time steps (or 10 fs), as in our previous work Q]. For REUS simulations, on the other hand, 
replica exchange was tried every 400 steps (or 200 fs), which is less frequent than in REMl. This is because we wanted 
to ensure sufficient time for system relaxation after A-parameter exchange. 

In REMl there are 16 replicas with 16 different temperatures listed in Table I. The temperatures are distributed 
exponentially, following the optimal distribution found in Ref. |4J|. After every 10 fs of parallel MD simulations, 
eight pairs of replicas corresponding to neighboring temperatures were simultaneously exchanged, and the pairing was 
alternated between the two possible choices Mm. 



For umbrella potentials, the 01 to H5 hydrogen-bonding distance, or "end-to-end distance," was chosen as the 
reaction coordinate ^ and the harmonic restraining potentials of ^ in Eq. ( p5| ) were imposed. The force constants, k^, 
and the midpoint positions, dg, are listed in Table I. 

In REUSl and USl, 14 replicas were simulated with the same set of umbrella potentials at T = 300 K. The first 
parameter value, 0.0 (0.0), in Table I means that the restraining potential is null, i.e., Vi = 0. The remaining 13 sets 
of parameters are the same as those adopted in Ref. [£l|. Let us order the umbrella potentials, Ve in Eq. (p4[), in the 
increasing order of the midpoint value di, i.e., the same order that appears in Table I. We prepared replicas so that 
the potential energy for each replica includes exactly one umbrella potential (here, we have M = L = 14). Namely, 
in Eq. (e3) for A = A^ we set 

A^^ = ^£,™ , (30) 

where dk^i is Kronecker's delta function, and we have 

E^Jq^^) = Eo{q^''^) + V^{q^^). (31) 

The difference between REUSl and USl is whether replica exchange is performed or not during the parallel MD 
simulations. In REUSl seven pairs of replicas corresponding to "neighboring" umbrella potentials, Vm and Vm+i, were 
simultaneously exchanged after every 200 fs of parallel MD simulations, and the pairing was alternated between the 
two possible choices. (Other pairings will have much smaller acceptance ratios of replica exchange.) The acceptance 
criterion for replica exchange is given by Eq. dlj), where Eq. (Hq) now reads (with the fixed inverse temperature 
P = l/SOOfcs) 

A = /? (k, (gl-'l) - V„, (gW) - V„,+i (gl-'l) + V,n+i (gl'l) ) , (32) 

where replica i and j respectively have umbrella potentials Vm and Vm+i before the exchange. 

In REUS2 and US2, 16 replicas were simulated at four different temperatures with four different restraining 
potentials (there are L — 4 umbrella potentials at Nt = 4 temperatures, making the total number of replicas 
M = Nt X i = 16; see Table I). We can introduce the following labeling for the parameters characterizing the 
replicas: 

Am = (Tm,Am) > Aj ,j — {Tj , Xj) . , ^ 

(m = l,-..,M) (/ = 1,-..,7Vt, J=l,---,i) ^^ 

The potential energy is given by Eq. (|3l|) with the replacement: m ^ J. Let us again order the umbrella potentials, 
Vj, and the temperatures, Tj, in the same order that appear in Table I. The difference between REUS2 and US2 
is whether replica exchange is performed or not during the MD simulations. In REUS2 we performed the following 
replica-exchange processes alternately after every 200 fs of parallel MD simulations: 

1. Exchange pairs of replicas corresponding to neighboring temperatures, T/ and Tf+i (i.e., exchange replicas i 
and j that respectively correspond to parameters Ajj and A/+i^j). (We refer to this process as T-exchange.) 

2. Exchange pairs of replicas corresponding to "neighboring" umbrella potentials, Vj and Vj+i (i.e., exchange 
replicas i and j that respectively correspond to parameters A/^j and A/.j+i). (We refer to this process as 
A- exchange.) 

In each of the above processes, two pairs of replicas were simultaneously exchanged, and the pairing was further 
alternated between the two possibilities. The acceptance criterion for these replica exchanges is given by Eq. ^u\), 
where Eq. (O) now reads 

A^iPi- (3i+,) (e^ (ql-'l) + Vj (gl^l) - E^ (gW) - Vj (gW)) , (34) 

for T-exchange, and 

A = /?/ (Vj (ql^'l) - Vj (gW) - Kz+i (g'-'O + Kz+i (g'^O) ' (35) 



for A-exchange. By this procedure, the random walk in the reaction coordinate space as well as in temperature space 
can be realized. Note that we carry out the velocity rescaling of Eq. (0) in T-exchange. In principle, we can also 
introduce a similar velocity rescaling in A-exchange to the two relevant atoms 01 and H5 in order to adjust for the 



exchange of the restraining potentials (because the restraining force acts only on 01 and H5). We also incorporated 
this rescaling but did not see much improvement in performance. The results presented below are thus those from no 
velocity rescaling in A-exchange. 

We now give the details of the results obtained in the present work. First of all, we examine whether the replica- 
exchange processes properly occurred in REM and REUS simulations. One criterion for the optimal performance is: 
whether the acceptance ratio of replica exchange is uniform and sufficiently large or not. In Tables II-IV we list the 
acceptance ratios of replica exchange corresponding to the adjacent pairs of temperatures or the restraining potentials. 
In all cases the acceptance ratios are almost uniform and large enough (> 10 %); all simulations indeed performed 
properly. In particular, the acceptance ratios for exchanging adjacent temperatures are significantly uniform in all 
cases, implying that the exponential temperature distributions of Ref. |Q are again optimal. However, the acceptance 
ratios for exchanging adjacent restraining potentials are not perfectly uniform, and there is some room for fine tuning. 
The acceptance ratio for exchanging restraining potentials depends on the strength of the force constants, ki, (see 
Eq. (p5|)) and we weakened the value from 1.2 kcal/mol-A^ in REUSl to 0.3 kcal/mol-A^ in REUS2 in order to have 
sufficient replica exchanges in REUS2. This is because we have a much smaller number of restraining potentials in 
REUS2 than in REUSl (3 versus 13), and yet both simulations have to cover the same range of reaction coordinate 
$, i.e., from A to 13.8 A. 

In order to have sufficient replica exchanges between neighboring temperatures and between neighboring restraining 
potentials, the probability distributions corresponding to neighboring parameters should have enough overlaps. In 
Fig. 1 (a) the canonical probability distributions of the unbiased potential energy Eq at the four chosen temperatures 
are shown. The results are for the parameters A/^i {I — 1, • • • ,4), i.e., for the case of no restraining potentials, and 
were obtained from the REUS2 simulation. In Fig. 1(b) the probability distributions of the reaction coordinate ^ with 
the four chosen restraining potentials are shown. The results are for the parameters A2^j {J — li ■ ■ ■ j4), i.e., for the 
temperature T = 315 K, and were also obtained from the REUS2 simulation. In both figures we do observe sufficient 
overlaps in pairs of the distributions corresponding to the neighboring parameter values, and this is reflected in the 
reasonable acceptance ratios listed in Table IV. 

In order to further confirm that our REM simulations performed properly, we have to examine the time series of 
various quantities and observe random walks. For instance, in Fig. 2 the trajectories of a few quantities in REUS2 are 
shown. In Fig. 2(a) we show the time series of replica exchange for the parameter Aii = (Ti, Ai) (i.e., Ti = 250 K 
and ki = di = 0.0). We do observe a random walk in replica space, and we see that all the replicas frequently visited 
the parameter value Ai j. 

The complementary picture to this is the time series of T-exchange and A-exchange for each replica. Free random 
walks both in "temperature space" and in "restraining potential space" were indeed observed. For instance, the time 
series of temperature exchange for one of the replicas (replica 1) is shown in Fig. 2(b). The corresponding time series 
of the reaction coordinate £,, the distance between atoms 01 and 115, for the same replica is shown in Fig. 2(c). We 
see that the conformational sampling along the reaction coordinate is significantly enhanced. In the blocked alanine 
trimer, the reaction coordinate ^ can be classified into three regions [|l|: the helical region (^ < 3 A), the turn region 
(3 A < ^ < 7 A), and the extended region (^ > 7 A). Thus, Fig. 2(c) implies that helix-coil transitions frequently 
occurred during the replica-exchange simulation, whereas in the conventional canonical simulations such a frequent 
folding and unfolding process cannot be seen. 

After confirming that the present REM and REUS simulations performed properly, we now present and compare 
the physical quantities calculated by these simulations. In Fig. 3 the potentials of mean force (PMF) of the unbiased 
system along the reaction coordinate ^ at T = 300 K are shown. The results are from REMl, REUSl, and USl 
simulations. For these calculations, the WHAM equations of Eqs. (|2^) and ( P7|) were solved by iteration first, and 
then Eq. ( p9| ) was used to obtain the PMF. We remark that for biomolecular systems the results obtained from the 
WHAM equations are insensitive to the values of gm in Eq. (g6|) |j]. Hence, we set gm = constant in the present 
article. From Fig. 3 we see that the PMF curves obtained by REMl and REUSl are essentially identical for low 
values of ^ (^ < 7 A). The two PMF curves start deviating slightly, as ^ gets larger, and for ^ > 9 A the agreement 
completely deteriorates. The disagreement comes from the facts that the average £, at the highest temperature in 
REMl (Tie = 1500 K) is < ^ >Tie— ^-O A and that the original REM with T-exchange only cannot sample accurately 
the region where ^ is much larger than < ^ >Tie ■ These two simulations were performed under very different conditions: 
One was run at different temperatures without restraining potentials and the other at one temperature with many 
restraining potentials (see Table I). We thus consider the results to be quite reliable for (^ < 9 A). 

On the other hand, the PMF obtained by USl is relatively larger than those obtained by REMl and REUSl in 
the region of 2 A < ^ < 4 A, which corresponds to the structural transition state between the a-helical and turn 
structures. This suggests that USl got trapped in states of energy local minima at T = 300 K. In the region of 
completely extended structures (^ > 9 A), the results of REUSl and USl are similar but the discrepancy is again 
non-negligible. We remark that at T = 300 K the PMF is the lowest for i^ = 2 A, which implies that the a-helical 
structure is favored at this temperature. 



We next study the temperature dependence of physical quantities obtained from the REMl, REUS2, and US2 
simulations. In Fig. 4(a) we show the PMF again at T = 300 K. We observe that the PMF curves from REMl and 
REUS2 are essentially identical for ^ < 9 A and that they deviate for ^ > 9 A, because the results for REMl are not 
reliable in this region as noted above. In fact, by comparing Figs. 3 and 4(a), we find that the PMF obtained from 
REUSl and REUS2 are almost in complete agreement dX T — 300 K in the entire range of ^ values shown. On the 
other hand, we observe a discrepancy between REUS2 and US2 results. The PMF curve for US2 is significantly less 
than that for REUS2 in the region 2 A < ^ < 8 A. Note that the PMF curves for USl and US2 are completely in 
disagreement (compare Figs. 3 and 4(a)). 

In Fig. 4(b) we show the PMF at T = 500 K, which we obtained from REMl, REUS2, and US2 simulations. We 
again observe that the results from REMl and REUS2 are in good agreement for a wide range of ^ values. We find 
that the results from REMl do not significantly deteriorate until ^> llAatT = 500 K, whereas it did start deviating 
badly for ^ > 9 A at T = 300 K. The PMF curve for US2 deviates strongly from the REUS2 results for ^ > 6 A and 
is much larger than that of REUS2 (and REMl) in this region. We remark that at T = 500 K the PMF is the lowest 
for ^ « 6 A, which implies that extended structures are favored at this temperature. 

In Fig. 5 we show the average values of the reaction coordinate ^ as a function of temperature. The results are again 
from the REMl, REUS2, and US2 simulations. The expectation values were calculated from Eq. (Eq). We find that 
the average reaction coordinate, or the average end-to-end distance, grows as the temperature is raised, reflecting the 
unfolding of the peptide upon increased thermal fluctuations. Again we observe an agreement between REMl and 
REUS2, whereas the results of US2 deviate. 

Let us emphasize that the total length of the MD simulations was the same (2 ns) for each replica in all the 
simulations performed. Hence, we have shown that the replica-exchange umbrella sampling can give much more 
accurate free-energy profiles along a reaction coordinate than the conventional umbrella sampling. 

IV. CONCLUSIONS 

In this article we have presented a multidimensional extension of the original replica-exchange method. One example 
of this approach is the combination of the replica-exchange method with the umbrella sampling, which we refer to 
as the replica- exchange umbrella sampling (REUS). While pairs of replicas with different temperatures are exchanged 
during the simulation in the original replica-exchange method, pairs of replicas with different temperatures and/or 
different biasing potentials for the umbrella sampling are exchanged in REUS. This greatly enhances the sampling of 
the conformational space and allows accurate calculations of free energy in a wide temperature range from a single 
simulation run, using the weighted histogram analysis method. The difference between REUS and the conventional 
umbrella sampling is just whether the replica-exchange process is performed or not. Only minor modifications to the 
conventional umbrella sampling method are necessary. However, the advantage of REUS over the umbrella sampling 
is significant, and the effectiveness was established with the system of an alanine trimer. 
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TABLE I. Summary of the replica parameters for the present simulations. 



Run" 



M" 



Nt 



Temperature [K] 



de [A] {ke [kcal/mol-A^])" 



REMl 



REUSl, USl 



REUS2, US2 



16 



14 



16 



16 



200, 229, 262, 299, 
342, 391, 448, 512, 
586, 670, 766, 876, 
1002, 1147, 1311, 
1500 
300 



250, 315, 397, 500 



14 0.0 (0.0)'', 1.8 (1.2), 2.8 (1.2), 3.8 (1.2), 

4.8 (1.2), 5.8 (1.2), 6.8 (1.2), 7.8 (1.2), 
8.8 (1.2), 9.8 (1.2), 10.8 (1.2), 11.8 (1.2), 
12.8 (1.2), 13.8 (1.2) 
4 0.0 (0.0), 7.8 (0.3), 10.8 (0.3), 13.8 (0.3) 



° REM, REUS, and US stand for an original replica-exchange simulation, replica-exchange umbrella sampling simu- 
lation, and conventional umbrella sampling simulation, respectively. 

^ M, Nt, and L are the total numbers of replicas, temperatures, and restraining potentials, respectively (see Eqs. (|lq ) 
and (|2J)). In REUS2 and US2 we set M = Nt x L for simplicity. We remark that this relation is not always 
required. For instance, the 16 replicas could have 16 different temperatures with 16 different restraining potentials 
(i.e., M = Nt = L = 16). 

■^ de and k^ {£ = 1, ■ • • , L) are the strengths and the midpoints of the restraining potentials, respectively (see Eq. (p5D 
'^ The parameter value 0.0 (0.0) means that the restraining potential is null, i.e., V^ — 0. 



TABLE II. Acceptance ratios of replica exchange in REMl. 



Pair of Temperatures 



Acceptance Ratio 



200 
229 
262 
299 
342 
391 
448 
512 
586 
670 
766 
876 ^ 
1002 
1147 
1311 



y 229 

> 262 

> 299 

> 342 

> 391 
>448 

> 512 

> 586 

> 670 

> 766 

> 876 
1002 

> 1147 

> 1311 

> 1500 



0.430 
0.433 
0.433 
0.428 
0.430 
0.423 
0.429 
0.427 
0.434 
0.437 
0.445 
0.446 
0.446 
0.454 
0.452 



11 



TABLE III. Acceptance ratios of replica exchange in REUSl. 



Pair of Restraint Parameters 



Acceptance Ratio 



0.0 (0.0) 


< — > 1.8 (1.2) 


0.202 


1.8 (1.2) 


< > 2.8 (1.2) 


0.210 


2.8 (1.2) 


< — > 3.8 (1.2) 


0.174 


3.8 (1.2) 


< > 4.8 (1.2) 


0.161 


4.8 (1.2) 


< — > 5.8 (1.2) 


0.223 


5.8 (1.2) 


< — > 6.8 (1.2) 


0.155 


6.8 (1.2) 


< — > 7.8 (1.2) 


0.211 


7.8 (1.2) 


< — > 8.8 (1.2) 


0.229 


8.8 (1.2) 


< — y 9.8 (1.2) 


0.119 


9.8 (1.2) 


< — > 10.8 (1.2) 


0.276 


10.8 (1.2) 


< — > 11.8 (1.2) 


0.156 


11.8 (1.2) 


< > 12.8 (1.2) 


0.138 


12.8 (1.2) 


< — > 13.8 (1.2) 


0.383 



TABLE IV. Acceptance ratios of replica exchange in REUS2. 



Temperature 



Pair of Restraint Parameters 



Acceptance Ratio 



250 
250 
250 
315 
315 
315 
397 
397 
397 
500 
500 
500 



0.0 (0.0) 

7.8 (0.3) < 

10.8 (0.3) 

0.0 (0.0) 

7.8 (0.3) < 

10.8 (0.3) 

0.0 (0.0) 

7.8 (0.3) < 

10.8 (0.3) 

0.0 (0.0) 

7.8 (0.3) < 

10.8 (0.3) 



► 7.8 (0.3) 
10.8 (0.3) 

► 13.8 (0.3) 

► 7.8 (0.3) 
10.8 (0.3) 

► 13.8 (0.3) 

► 7.8 (0.3) 
10.8 (0.3) 

► 13.8 (0.3) 

► 7.8 (0.3) 
10.8 (0.3) 

► 13.8 (0.3) 



0.149 
0.104 
0.127 
0.250 
0.105 
0.120 
0.363 
0.135 
0.160 
0.483 
0.185 
0.228 



Restraint Parameters 



Pair of Temperatures 



Acceptance Ratio 



0.0 (0.0) 
0.0 (0.0) 
0.0 (0.0) 
7.8 (0.3) 
7.8 (0.3) 
7.8 (0.3) 
10.8 (0.3) 
10.8 (0.3) 
10.8 (0.3) 
13.8 (0.3) 
13.8 (0.3) 
13.8 (0.3) 



250 < — » 315 
315 < — > 397 
397 < — > 500 
250 < — > 315 
315 < — » 397 
397 < — > 500 
250 < — > 315 
315 < — > 397 
397 < — > 500 
250 < — > 315 
315 < — > 397 
397 < — > 500 



0.193 
0.186 
0.189 
0.174 
0.179 
0.190 
0.189 
0.184 
0.195 
0.185 
0.184 
0.205 
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Figure Captions 

• Fig. 1. Probability distributions obtained from the replica-exchange umbrella sampling simulation (REUS2). 
(a) The canonical probability distributions of the unbiased potential energy Eq at the four chosen temperatures 
(the curves from left to right correspond to T = 250, 315, 397, 500 K). The results are for the parameters A/^i 
(/ = 1, • • • , 4), i.e., for the case of no restraining potentials (see Table I), (b) The probability distributions of the 
reaction coordinate ^, the distance between the atoms 01 and H5, with the four chosen restraining potentials 
(the curves from left to right correspond to de [A] {ki [kcal/mol-A^]) = 0.0 (0.0), 7.8 (0.3), 10.8 (0.3), 13.8 (0.3)). 
The results are for the parameters A2.,/ {J ~ 1, • • • ,4), i.e., for the temperature T = 315 K (see Table I). 

• Fig. 2. Time series of the replica-exchange umbrella sampling simulation (REUS2). (a) Replica exchange for 
the parameter Ai_i = (Ti,Ai) (i.e., Ti = 250 K and ki = di = 0.0, see Table I), (b) Temperature exchange for 
one of the replicas (replica 1). (c) The reaction coordinate ^ for one of the replicas (replica 1). 

• Fig. 3. The PMF along the reaction coordinate ^ at T = 300 K. The dotted, solid, and dashed curves were ob- 
tained from the original REM (REMl), the replica-exchange umbrella sampling (REUSl), and the conventional 
umbrella sampling (USl), respectively. 

• Fig. 4. The PMF along the reaction coordinate £, at two temperatures, (a) The PMF at T = 300 K. The dotted, 
solid, and dashed curves were obtained from the original REM (REMl), the replica-exchange umbrella sampling 
(REUS2), and the conventional umbrella sampling (US2), respectively, (b) The PMF at T = 500 K. The dotted, 
solid, and dashed curves were obtained from the original REM (REMl), the replica-exchange umbrella sampling 
(REUS2), and the conventional umbrella sampling (US2), respectively. 

• Fig. 5. Average values of the reaction coordinate ^ as a function of temperature. The dotted, solid, and dashed 
curves were obtained from the original REM (REMl), the replica-exchange umbrella sampling (REUS2), and 
the conventional umbrella sampling (US2), respectively. Although the highest temperature in REMl is 1500 K, 
only the results for the temperature range between 200 K and 1000 K are shown for REMl. Since the lowest 
and highest tempeatures in REUS2 and US2 are respectively 250 K and 500 K, only the results between these 
temperatures are shown for these simulations. 
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